Antje Hofmann: Learning Syntactic Relations with the Mole Task

RePsychLing in SMLP2022

1 Background

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

• Original analysis is by Antje Hofmann.
• MixedModels.jl version
• Addition of new chunks illustrate
  • selection of parsimonious LMM using random-effects PCA
  • plotting conditional means
  • illustration of borrowing strength

2 Readme

2.1 Design

2.2 Variables

• Subj: Participant ID (renamed from ID; random factor)
• Item: Word ID (random factor)
• age: 4 - 8 years
• Block (within-Subj/within-Item):
  • 1st Learning
  • 2nd Learning
  • Disruption
  • Recovery
• Target(renamend fom targetness)
  • non-syllable target
  • syllable target
• rt: response time

3 Setup

3.1 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code

using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");

• The data are available as an arrow file.
• Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
• Order of factor levels should be checked.

Code

dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat,
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age (linear) at 6 years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)

8 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	Type
1	Subj		p01		p85	0	Union{Missing, String}
2	Item		01PA.ogg		32SG.ogg	0	Union{Missing, String}
3	age	6.21081	3.7	6.31	8.04	0	Union{Missing, Float64}
4	Block		1st Learning		Recovery	0	Union{Missing, CategoricalValue{String, UInt32}}
5	Target		Non-target syllable		Target syllable	0	Union{Missing, CategoricalValue{String, UInt32}}
6	rt	2335.39	228.0	2100.0	8389.0	0	Union{Missing, Float64}
7	a1	0.210813	-2.3	0.31	2.04	0	Float64
8	lrt	7.67367	5.42935	7.64969	9.03468	0	Float64

• Centering age at six years yields an interpretable GM
• Factor levels can also be set when contrasts are defined (see below).
• BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
• Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code

contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code

f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)

Minimizing 65    Time: 0:00:01 (17.57 ms/it)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.6895	0.0371	207.05	<1e-99	0.0239	0.2637
Block: 2nd Learning	-0.0920	0.0104	-8.87	<1e-18
Block: Disruption	0.0516	0.0126	4.09	<1e-04
Block: Recovery	-0.0477	0.0125	-3.81	0.0001
Target: Target syllable	-0.0108	0.0040	-2.69	0.0072
a1	-0.0896	0.0324	-2.77	0.0056
Block: 2nd Learning & Target: Target syllable	-0.0057	0.0104	-0.55	0.5854
Block: Disruption & Target: Target syllable	-0.0262	0.0124	-2.11	0.0347
Block: Recovery & Target: Target syllable	0.0020	0.0123	0.16	0.8741
Block: 2nd Learning & a1	0.0347	0.0091	3.83	0.0001
Block: Disruption & a1	0.0214	0.0108	1.97	0.0486
Block: Recovery & a1	-0.0159	0.0108	-1.47	0.1409
Target: Target syllable & a1	-0.0052	0.0035	-1.46	0.1430
Block: 2nd Learning & Target: Target syllable & a1	-0.0113	0.0090	-1.25	0.2111
Block: Disruption & Target: Target syllable & a1	0.0113	0.0108	1.05	0.2933
Block: Recovery & Target: Target syllable & a1	0.0086	0.0107	0.80	0.4248
Residual	0.3004

4.3 Extract indicator variables

Code

X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	Type
1	Subj		p01		p85	0	Union{Missing, String}
2	Item		01PA.ogg		32SG.ogg	0	Union{Missing, String}
3	age	6.21081	3.7	6.31	8.04	0	Union{Missing, Float64}
4	Block		1st Learning		Recovery	0	Union{Missing, CategoricalValue{String, UInt32}}
5	Target		Non-target syllable		Target syllable	0	Union{Missing, CategoricalValue{String, UInt32}}
6	rt	2335.39	228.0	2100.0	8389.0	0	Union{Missing, Float64}
7	a1	0.210813	-2.3	0.31	2.04	0	Float64
8	lrt	7.67367	5.42935	7.64969	9.03468	0	Float64
9	trng	-0.014562	-0.75	0.25	0.25	0	Float64
10	drpt	-0.0528019	-0.5	-0.5	0.5	0	Float64
11	rcvr	0.0524467	-0.25	-0.25	0.75	0	Float64
12	trgt	0.00331492	-1.0	1.0	1.0	0	Float64

Switch to indicator variables and refit the model.

Code

f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.6895	0.0371	207.05	<1e-99	0.0239	0.2637
trng	-0.0920	0.0104	-8.87	<1e-18
drpt	0.0516	0.0126	4.09	<1e-04
rcvr	-0.0477	0.0125	-3.81	0.0001
trgt	-0.0108	0.0040	-2.69	0.0072
a1	-0.0896	0.0324	-2.77	0.0056
trng & trgt	-0.0057	0.0104	-0.55	0.5854
drpt & trgt	-0.0262	0.0124	-2.11	0.0347
rcvr & trgt	0.0020	0.0123	0.16	0.8741
trng & a1	0.0347	0.0091	3.83	0.0001
drpt & a1	0.0214	0.0108	1.97	0.0486
rcvr & a1	-0.0159	0.0108	-1.47	0.1409
trgt & a1	-0.0052	0.0035	-1.46	0.1430
trng & trgt & a1	-0.0113	0.0090	-1.25	0.2111
drpt & trgt & a1	0.0113	0.0108	1.05	0.2933
rcvr & trgt & a1	0.0086	0.0107	0.80	0.4248
Residual	0.3004

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code

f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)

Minimizing 463   Time: 0:00:00 ( 1.23 ms/it)

false

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0742252	0.2724431
	Block: 2nd Learning	0.0362305	0.1903431	.
	Block: Disruption	0.0156482	0.1250926	.	.
	Block: Recovery	0.0133091	0.1153650	.	.	.
	Target: Target syllable	0.0002450	0.0156523	.	.	.	.
	Block: 2nd Learning & Target: Target syllable	0.0001921	0.0138594	.	.	.	.	.
	Block: Disruption & Target: Target syllable	0.0032728	0.0572082	.	.	.	.	.	.
	Block: Recovery & Target: Target syllable	0.0006683	0.0258522	.	.	.	.	.	.	.
Item	(Intercept)	0.0009504	0.0308293
	a1	0.0006760	0.0259993	-0.83
Residual		0.0767202	0.2769841

Again, check the equivalence.

Code

f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  +
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)

Minimizing 463   Time: 0:00:00 ( 0.53 ms/it)

false

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0742252	0.2724431
	trng	0.0362305	0.1903431	.
	drpt	0.0156482	0.1250926	.	.
	rcvr	0.0133091	0.1153650	.	.	.
	trgt	0.0002450	0.0156523	.	.	.	.
	trng & trgt	0.0001921	0.0138594	.	.	.	.	.
	drpt & trgt	0.0032728	0.0572082	.	.	.	.	.	.
	rcvr & trgt	0.0006683	0.0258522	.	.	.	.	.	.	.
Item	(Intercept)	0.0009504	0.0308293
	a1	0.0006760	0.0259993	-0.83
Residual		0.0767202	0.2769841

4.5 A complex parameter LMM m_cpx

Code

m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx))
VarCorr(m_cpx)

Minimizing 1355      Time: 0:00:00 ( 0.64 ms/it)

true

Principal components based on correlation matrix

 (Intercept)   1.0     .      .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .      .
 trng          0.35  -0.35   1.0     .      .      .      .      .
 drpt          0.06  -0.06   0.35   1.0     .      .      .      .
 rcvr         -0.42   0.42  -0.23  -0.48   1.0     .      .      .
 trgt & trng   0.26  -0.26  -0.12   0.09  -0.16   1.0     .      .
 trgt & drpt  -0.4    0.4   -0.07   0.14  -0.01  -0.76   1.0     .
 trgt & rcvr   0.61  -0.61   0.16   0.04  -0.3    0.48  -0.85   1.0

Normalized cumulative variances:

[0.441, 0.667, 0.8072, 0.906, 0.9581, 1.0, 1.0, 1.0]

Component loadings
 

               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.46   0.14   0.36  -0.2   -0.28   0.12   0.05  -0.71
 trgt          0.46  -0.14  -0.36   0.2    0.28  -0.12  -0.05  -0.71
 trng         -0.19   0.43   0.17   0.76  -0.0   -0.4   -0.13   0.0
 drpt         -0.09   0.49  -0.59   0.15  -0.19   0.57   0.13  -0.0
 rcvr          0.26  -0.42   0.31   0.47  -0.39   0.51  -0.15   0.0
 trgt & trng  -0.31  -0.36  -0.49  -0.01  -0.52  -0.34  -0.38   0.0
 trgt & drpt   0.39   0.44   0.14  -0.29  -0.16   0.02  -0.72   0.0
 trgt & rcvr  -0.46  -0.18  -0.04   0.08   0.6    0.33  -0.53   0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2393.5710                         
[2]         56  2353.2482  40.3229      28  0.0619
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0757920	0.2753034
	trgt	0.0001016	0.0100800	-1.00
	trng	0.0334990	0.1830274	+0.35	-0.35
	drpt	0.0198176	0.1407751	+0.06	-0.06	+0.35
	rcvr	0.0191682	0.1384494	-0.42	+0.42	-0.23	-0.48
	trgt & trng	0.0020614	0.0454021	+0.26	-0.26	-0.12	+0.09	-0.16
	trgt & drpt	0.0145492	0.1206203	-0.40	+0.40	-0.07	+0.14	-0.01	-0.76
	trgt & rcvr	0.0085469	0.0924494	+0.61	-0.61	+0.16	+0.04	-0.30	+0.48	-0.85
Item	(Intercept)	0.0000000	0.0000000
	a1	0.0004752	0.0217991	+NaN
Residual		0.0768486	0.2772159

The deviance improves, but we end up with an overparameterized LMM. Remove the two small VCs of the three Block x Target interactions (i.e., the third contrast of the Block factor).

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code

m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # not ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx))
VarCorr(m_prm1)

Minimizing 1465      Time: 0:00:00 ( 0.61 ms/it)

true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -0.58   1.0     .      .      .      .      .
 drpt          0.04  -0.23   1.0     .      .      .      .
 rcvr         -0.41   0.21  -0.47   1.0     .      .      .
 trgt & drpt  -0.42   0.3    0.21  -0.06   1.0     .      .
 trgt & rcvr   0.65  -0.22   0.03  -0.29  -0.87   1.0     .
 trng          0.33  -0.32   0.32  -0.21  -0.12   0.22   1.0

Normalized cumulative variances:
[0.4076, 0.6504, 0.7761, 0.8828, 0.9634, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.5    0.02  -0.05  -0.36  -0.38  -0.62  -0.3
 trgt          0.38   0.14   0.56   0.49  -0.35  -0.3   -0.26
 drpt         -0.15  -0.61   0.15   0.25   0.61  -0.36  -0.14
 rcvr          0.29   0.44  -0.56   0.23   0.27  -0.52   0.1
 trgt & drpt   0.4   -0.49  -0.05  -0.23  -0.33  -0.3    0.59
 trgt & rcvr  -0.49   0.29   0.35   0.24   0.07  -0.13   0.68
 trng         -0.3   -0.29  -0.47   0.64  -0.42   0.15  -0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         48  2334.2053   59.3657      20  <1e-05
[3]         56  2353.2482  -19.0428       8     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0759692	0.2756251
	trgt	0.0003130	0.0176929	-0.58
	drpt	0.0194377	0.1394192	+0.04	-0.23
	rcvr	0.0188254	0.1372057	-0.41	+0.21	-0.47
	trgt & drpt	0.0109856	0.1048121	-0.42	+0.30	+0.21	-0.06
	trgt & rcvr	0.0075846	0.0870893	+0.65	-0.22	+0.03	-0.29	-0.87
	trng	0.0354738	0.1883450	+0.33	-0.32	+0.32	-0.21	-0.12	+0.22
Item	(Intercept)	0.0008077	0.0284198
	a1	0.0006179	0.0248571	-0.94
Residual		0.0761744	0.2759972

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code

m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx))
VarCorr(m_prm2)

Minimizing 740   Time: 0:00:00 ( 0.55 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.22   1.0     .      .     .
 drpt & trgt  -0.36   0.26  -0.1   -0.07   1.0     .     .
 rcvr & trgt   0.59  -0.0   -0.27   0.12  -0.87   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3444, 0.5975, 0.7404, 0.8569, 0.9478, 0.9939, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.07  0.0   -0.04   0.61  -0.58   0.15
 drpt         -0.09   0.61  0.0   -0.05  -0.62  -0.47   0.11
 rcvr          0.3   -0.49  0.0    0.52  -0.16  -0.58  -0.17
 trng         -0.26   0.36  0.0    0.83   0.09   0.3   -0.08
 drpt & trgt   0.47   0.45  0.0   -0.1    0.37  -0.14  -0.65
 rcvr & trgt  -0.59  -0.22  0.0   -0.14  -0.27   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         42  2341.3452   52.2258      14  <1e-05
[3]         56  2353.2482  -11.9030      14     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0760175	0.2757128
	drpt	0.0193591	0.1391371	+0.04
	rcvr	0.0188566	0.1373195	-0.41	-0.47
	trng	0.0351279	0.1874243	+0.33	+0.34	-0.22
	drpt & trgt	0.0099356	0.0996773	-0.36	+0.26	-0.10	-0.07
	rcvr & trgt	0.0066342	0.0814504	+0.59	-0.00	-0.27	+0.12	-0.87
	trgt	0.0002660	0.0163107	.	.	.	.	.	.
Item	(Intercept)	0.0008269	0.0287560
	a1	0.0006738	0.0259576	-0.93
Residual		0.0762049	0.2760523

Check the Item-related CP. It is very large, might be spurious.

Code

m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx))
VarCorr(m_prm3)

Minimizing 1162      Time: 0:00:00 ( 0.53 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.21   1.0     .      .     .
 drpt & trgt  -0.37   0.25  -0.08  -0.07   1.0     .     .
 rcvr & trgt   0.61   0.04  -0.31   0.13  -0.86   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3499, 0.5979, 0.7408, 0.8584, 0.9488, 0.9943, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.03  0.0   -0.01   0.61  -0.57   0.16
 drpt         -0.12   0.62  0.0   -0.06  -0.61  -0.47   0.12
 rcvr          0.33  -0.47  0.0    0.52  -0.19  -0.58  -0.18
 trng         -0.27   0.35  0.0    0.83   0.08   0.31  -0.08
 drpt & trgt   0.45   0.47  0.0   -0.1    0.38  -0.15  -0.63
 rcvr & trgt  -0.58  -0.22  0.0   -0.15  -0.26   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         41  2356.2006   37.3704      13  0.0004
[3]         42  2341.3452   14.8554       1  0.0001
[4]         56  2353.2482  -11.9030      14     NaN
───────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0760533	0.2757776
	drpt	0.0193959	0.1392690	+0.04
	rcvr	0.0189445	0.1376392	-0.41	-0.47
	trng	0.0350939	0.1873337	+0.33	+0.34	-0.21
	drpt & trgt	0.0091654	0.0957362	-0.37	+0.25	-0.08	-0.07
	rcvr & trgt	0.0059703	0.0772677	+0.61	+0.04	-0.31	+0.13	-0.86
	trgt	0.0002813	0.0167707	.	.	.	.	.	.
Item	(Intercept)	0.0006679	0.0258444
	a1	0.0005355	0.0231408	.
Residual		0.0763232	0.2762666

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code

let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum ∘ leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

	model	pars	geomdof	AIC	AICc	BIC
	Symbol	Int64	Int64	Int64	Int64	Int64
1	m_zcp2	28	292	2450	2450	2639
2	m_prm2	42	287	2425	2426	2709
3	m_prm1	48	282	2430	2431	2754
4	m_cpx	56	259	2465	2466	2843

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code

coeftable(m_prm2)

	Coef.	Std. Error	z	Pr(>
(Intercept)	7.68255	0.0388114	197.95	<1e-99
trgt	-0.00948337	0.00442408	-2.14	0.0321
trng	-0.0841462	0.0279289	-3.01	0.0026
drpt	0.042879	0.022888	1.87	0.0610
rcvr	-0.0370272	0.0226406	-1.64	0.1020
a1	-0.0922378	0.0339009	-2.72	0.0065
trgt & trng	-0.00176196	0.00963308	-0.18	0.8549
trgt & drpt	-0.0199346	0.0181564	-1.10	0.2722
trgt & rcvr	-0.00411183	0.0162151	-0.25	0.7998
trgt & a1	-0.00874189	0.00389127	-2.25	0.0247
trng & a1	0.0414062	0.0244065	1.70	0.0898
drpt & a1	0.0287386	0.0198962	1.44	0.1486
rcvr & a1	-0.0208128	0.0197094	-1.06	0.2910
trgt & trng & a1	-0.00234866	0.00840927	-0.28	0.7800
trgt & drpt & a1	0.017156	0.0158645	1.08	0.2795
trgt & rcvr & a1	0.00435148	0.0141861	0.31	0.7590

4.8 Pruning fixed-effects

There is an option of pruning some higher-order interactions.

4.9 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.9.1 Residuals over fitted

Code

scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.9.2 Q-Q plot

Code

qqnorm(m_prm2; qqline=:none)

Hm?

4.9.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code

let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.10 Conditional means of random effects

4.10.1 Subject-related conditional means of random effects

Code

cm_Subj = first(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Subj; orderby=1)

Figure 2: Prediction intervals on subject random effects for model m_prm2

4.10.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code

shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

There is an outlier subject

Code

cm = raneftables(m_prm2);
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

52 rows × 8 columns (omitted printing of 1 columns)

	Subj	(Intercept)	drpt	rcvr	trng	drpt & trgt	rcvr & trgt
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	p10	-0.665239	-0.0800783	0.176854	-0.502966	0.3251	-0.281851
2	p14	-0.517129	0.149376	-0.104141	-0.380418	0.0865924	-0.0397509
3	p76	0.0792905	-0.100487	0.10478	-0.309599	-0.0284471	0.0204237
4	p02	-0.622029	-0.197088	0.260561	-0.251622	0.0703971	-0.126565
5	p36	0.421409	-0.0363361	0.0871325	-0.231898	-0.168933	0.145207
6	p39	0.291872	-0.130734	0.143494	-0.200954	-0.0192091	0.00241466
7	p74	-0.329558	-0.0609616	0.0915499	-0.179366	0.0670593	-0.0858263
8	p67	0.0383289	0.00117737	-0.0194406	-0.15504	0.0331733	-0.0118186
9	p06	-0.237713	0.139293	0.0263688	-0.117255	-0.0667732	0.0402579
10	p54	0.313212	-0.161867	0.0464652	-0.116672	-0.0367807	0.0410875
11	p23	-0.0288518	-0.0711039	0.0526805	-0.103828	-0.126974	0.0803409
12	p42	-0.0973798	0.288202	-0.16646	-0.103743	0.121444	-0.0338319
13	p59	0.421684	0.0211085	-0.255639	-0.102933	-0.0122904	0.0649662
14	p45	-0.363704	-0.295772	0.169706	-0.0975617	0.0434683	-0.0933397
15	p30	-0.32999	-0.150138	-0.0896487	-0.0878293	-0.0134365	-0.00936306
16	p83	-0.23045	-0.0533087	-0.160339	-0.0702177	0.0131955	-0.0083732
17	p84	-0.389921	-0.202987	0.444955	-0.0490962	-0.199401	0.0206631
18	p72	0.176323	-0.0684158	-0.0524196	-0.0387423	-0.119473	0.117996
19	p12	-0.148151	-0.0448463	0.0523027	-0.0284075	-0.0671666	0.0361057
20	p04	0.138257	-0.0447661	0.0390786	-0.027842	-0.0032299	0.00440696
21	p77	0.134281	-0.0964947	0.0522384	-0.021795	-0.0800806	0.0499684
22	p69	0.104395	-0.0390978	-0.0831703	-0.0144483	0.0226967	0.0349963
23	p79	-0.0769208	-0.0193177	0.0301171	-0.0144253	-0.0321678	0.0103847
24	p55	-0.198413	-0.0638211	-0.0784104	-0.00928064	-0.0624524	0.0241703
25	p27	0.0253274	-0.00351493	-0.00407116	-0.00610536	-0.0976157	0.0736163
26	p28	0.150447	-0.00600531	-0.151517	0.00701752	0.0683197	-0.0247804
27	p25	0.267637	0.137046	-0.0474272	0.00967245	-0.035619	0.0451562
28	p64	0.268335	0.00238078	0.0338587	0.0132533	-0.018809	0.024563
29	p57	0.00282995	-0.00253166	-0.0153407	0.0155963	0.0551781	-0.0361016
30	p53	-0.0783682	0.0257413	0.0203331	0.0180285	0.0430858	-0.0321112
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

It appears to be Subj “p10”.

4.10.3 Item-related conditional means of random effects

Code

cm_Item = last(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Item; orderby=1)

Figure 4: Prediction intervals on item random effects for model m_prm2

5 Appendix

Code

versioninfo()

Julia Version 1.8.1
Commit afb6c60d69a (2022-09-06 15:09 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 8 on 8 virtual cores